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Abstract. We investigate a prominent vibrational feature in amorphous silica, the 
so-called boson peak, by means of molecular dynamics computer simulations. The 
dynamic structure factor S(q, v) in the liquid, as well as in the glass state, scales roughly 
with temperature, in agreement with the harmonic approximation. By varying the size 
of the system and the masses of silicon and oxygen we show that the excitations giving 
rise to the boson peak are due to the coupling to transverse acoustic modes. 



I INTRODUCTION 

In the last few years various scattering techniques, such as neutron, Raman 
and X-ray scattering, have been used to investigate the so-called boson peak, a 
vibrational feature, which is found in the frequency spectra of many, typically 
strong, glass formers at a frequency of about 1 THz [1]. In this context various 
mechanisms giving rise to this peak have been proposed, such as certain localized 
vibrational modes or scattering of acoustic waves, and also simple models have 
been developed that produce an excess over the Debye behavior in the density of 
states [2]. 

Especially in the case of silica molecular dynamics computer simulations have 
recently been used in order to gain insight into the nature of the boson peak [3-5]. 
Despite the limitations of these simulations, such as the small system size (of the 
order of 10 3 -10 4 particles) and high cooling rates (of the order of 10 12 K/s), they 
are very useful because they include in principle the full microscopic information 
in form of the particle trajectories. Most of the recent computer simulation studies 
have investigated the boson peak within the harmonic approximation in that the 
eigenvalues and eigenvectors of the dynamical matrix have been calculated [3,4]. 
In contrast to this method we use the full microscopic information to determine 
quantities like the dynamic structure factor S(q, u) directly from the particle coor- 
dinates. Thus, we are not restricted to the harmonic approximation and we are able 



to compare the dynamics of our silica model in the liquid state with the dynamics 
in the glass state. Moreover, by varying parameters like the size of the system and 
the mass of the particles we gain information on the character of the boson peak 
excitations. 

II DETAILS OF THE SIMULATION 

The silica model we use for our simulation is the one proposed by van Beest et 
al. [6] which is given by 
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The values of the partial charges and the constants Aij, and can be found 
in the original publication. The simulations were done at constant volume keeping 
the density fixed at 2.37 g/cm 3 . Our simulation box contains 8016 particles with a 
box length of 48.37 A. We investigate the equilibrium dynamics of the liquid state 
as well as the glass state. The lowest temperature for which we were able to fully 
equilibrate our system was 2750 K. At this temperature we integrated the equations 
of motion over 13 million time steps of 1.6 fs, thus over a time span of about 21 ns. 
The glass state was produced by starting from two equilibrium configurations at 
T = 2900 K and cooling them to the temperatures T = 1670 K, 1050 K and 300 K 
with a cooling rate of 1.8 • 10 12 K/s. The details of how we calculated the time 
Fourier transformations can be found elsewhere [7]. 

Ill RESULTS 

We investigate the high frequency dynamics of silica by means of the dynamic 
structure factor 
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and its self part S s (q, v) which can be extracted from Eq. (2) by taking into account 
only the terms with k = I. In the following we will consider only S(q, v) for 
the oxygen-oxygen correlations because the oxygen-silicon and the silicon-silicon 
correlations behave similarly with respect to the features which are discussed below. 

As we have reported elsewhere for the liquid state at the temperature T = 
2900 K [5], apart from optical modes with frequencies v > 20 THZ, two types of 
excitations are visible in the dynamic structure factor for q > 0.23 A -1 . The first 
one corresponds to the boson peak which is located, essentially independent of q, 
around 1.8 THz. The second one corresponds to dispersive longitudinal acoustic 
modes. Note that the latter are not like longitudinal acoustic excitations in har- 
monic crystals because, due to the disorder, they cannot be described as plane 
waves. 
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FIGURE 1. S(q, v)/T as a function of frequency for different temperatures. 

Having found the aforementioned two features at T = 2900 K we will now look 
at the temperature dependence of S{q,v), which is shown in Fig. 1, by plotting 
S(q,u)/T versus v at q — 1.7 A -1 in the frequency range below 20 THz. From 
this figure we can conclude that the dynamic structure factor scales roughly with 
temperature which is expected if the harmonic approximation is valid. Moreover 
we can clearly identify for all temperatures two peaks: the boson peak located at 
1.75 THz (vertical line) and a peak which corresponds to the longitudinal acoustic 
modes located around 17 THz. Even at T = 3760 K the excitations giving rise to a 
boson peak at lower temperatures are at least partially present in that a shoulder 
can be recognized in the frequency region of the boson peak. 

In order to get some insight into the properties of the vibrational modes of our 
silica system at frequencies around 1 THz we varied the size of the system at a fixed 
mass density of 2.37 g/cm 3 . Fig. 2 shows the self part of the dynamic structure 
factor for N = 336, 1002 and 8016 particles at the temperature T = 3760 K 
and the three q values 0.37 A -1 , 1.7 A -1 and 4.75 A -1 . Whereas the curves for 
the different system sizes coincide for frequencies that are larger than a weakly iV 
dependent frequency Vcut(N), for v < ^cut(^) the amplitude of S s (q,u) decreases 
with decreasing N. Note that f cu t(-^0 is essentially independent of the wave-vector 
q. We read off z/ cut « 1.7 THz for N = 336 and z/ cut w 1.2 THz for N = 1002. Both 
frequencies are marked as vertical lines in Fig. 2. u cu t(N) coincides approximately 
with the frequency of the transverse acoustic excitation corresponding to the lowest 
q value which is determined by the size of the simulation box. To see that this is 
the case note that the lowest q values for N = 336 and N = 1002 are q min = 
0.37 A" 1 and q min = 0.26 A -1 , respectively. In comparison to that the q values 
we read off from the transverse acoustic dispersion branch for T = 3760 K and 
N = 8016, and which correspond to the frequency u cut for N = 336 and N = 1002, 
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FIGURE 2. Self part of the dynamic structure factor as a function of frequency for the system 
sizes N = 336, 1002 and 8016 at the temperature T = 3760 K. For the explanation of the vertical 
lines see text. 

are 0.32 A -1 and 0.22 A -1 , respectively (see [5]). That the latter q values are 
slightly smaller than the corresponding values for q min is due to the fact that the 
transverse dispersion branch has been determined from the peak maxima v max (q) 
in the transverse current correlation function. So there is always a significant 
contribution of transverse acoustic modes with frequencies v < f ma x(<?) for a given 
q. All this can be summarized by saying that the absence of transverse acoustic 
modes is connected with a missing of excitations giving rise to the boson peak. 
Therefore, in the smaller systems only the high frequency part of the boson peak is 
present. Moreover, it seems that the boson peak modes are only fully present, for a 
given frequency, if there exist transverse acoustic excitations at the same frequency. 

To learn more about the character of the boson peak excitations we varied also 
the masses of the silicon and oxygen atoms such that the mass density remains 
fixed. Fig. 3 shows S(q, v) at T — 2750 K and q = 0.6 A -1 for the four mass 
pairs Mi = (28.086,15.999), M 2 = (14.043,23.021), M 3 = (44.085,8.000), and 
M 4 = (56.085, 2.000) where the first and the second number are the masses in 
atomic units for silicon and oxygen, respectively. Note that M 1 corresponds to the 
real masses of silicon and oxygen normally used in our simulation. From the figure 
we see that there is a strong dependence on the mass ratio for the two peaks visible 
for Mi above 20 THz which are due to localized optical modes. In contrast to that, 
at least within the accuracy of the statistics of our data, there is no dependence 
for the modes giving rise to the boson peak and the acoustic modes which means 
that the boson peak excitations cannot be strongly localized. This supports the 
aforementioned statement that the boson peak is due to the coupling to transverse 
acoustic modes. The fact that the boson peak is independent of the mass ratio 
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FIGURE 3. Dynamic structure factor at T = 2750 K and q = 0.6 A 1 under variation of the 
masses of the silicon and oxygen atoms such that the mass density is fixed. 

between the silicon and oxygen mass is also a good test for theoretical models of 
the boson peak in silica. 

In conclusion we investigated the excitations giving rise to the boson peak by 
means of molecular dynamics computer simulations. We find that the dynamic 
structure factor S(q, v) scales roughly with temperature in the range 3760 K > 
T > 300 K, which means that our silica system is in this sense quite harmonic even 
for temperatures as high as 3760 K. By calculating S s (q,u) for different system 
sizes we find that the modes contributing to the boson peak are only fully present 
at a given frequency if there exist transverse acoustic modes at the same frequency. 
This is supported by the fact that the height and the width of the boson peak 
are independent under the variation of the masses of silicon and oxygen if the 
mass density is fixed. So we observe that the boson peak is due to a coupling to 
transverse acoustic modes. Of course, the explanation of the nature of this coupling 
is an interesting goal for the future. 
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